{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 13\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code from previous chapters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_system`, `plot_results`, and `calc_total_infected` are unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(beta, gamma):\n",
    "    \"\"\"Make a system object for the SIR model.\n",
    "    \n",
    "    beta: contact rate in days\n",
    "    gamma: recovery rate in days\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(S=89, I=1, R=0)\n",
    "    init /= np.sum(init)\n",
    "\n",
    "    t0 = 0\n",
    "    t_end = 7 * 14\n",
    "\n",
    "    return System(init=init, t0=t0, t_end=t_end,\n",
    "                  beta=beta, gamma=gamma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_results(S, I, R):\n",
    "    \"\"\"Plot the results of a SIR model.\n",
    "    \n",
    "    S: TimeSeries\n",
    "    I: TimeSeries\n",
    "    R: TimeSeries\n",
    "    \"\"\"\n",
    "    plot(S, '--', label='Susceptible')\n",
    "    plot(I, '-', label='Infected')\n",
    "    plot(R, ':', label='Recovered')\n",
    "    decorate(xlabel='Time (days)',\n",
    "             ylabel='Fraction of population')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_total_infected(results):\n",
    "    \"\"\"Fraction of population infected during the simulation.\n",
    "    \n",
    "    results: DataFrame with columns S, I, R\n",
    "    \n",
    "    returns: fraction of population\n",
    "    \"\"\"\n",
    "    return get_first_value(results.S) - get_last_value(results.S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Write a slope function for the SIR model and test it with `run_ode_solver`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def slope_func(state, t, system):\n",
    "    \"\"\"Update the SIR model.\n",
    "    \n",
    "    state: State (s, i, r)\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: State (sir)\n",
    "    \"\"\"\n",
    "    s, i, r = state\n",
    "\n",
    "    infected = system.beta * i * s    \n",
    "    recovered = system.gamma * i\n",
    "    \n",
    "    dSdt = -infected\n",
    "    dIdt = infected - recovered\n",
    "    dRdt = recovered\n",
    "    \n",
    "    return dSdt, dIdt, dRdt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.003658888888888889, 0.0008811111111111112, 0.002777777777777778)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = make_system(0.333, 0.25)\n",
    "slope_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sol</th>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_events</th>\n",
       "      <td>[]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nfev</th>\n",
       "      <td>212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>njev</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nlu</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>status</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "sol                                                      None\n",
       "t_events                                                   []\n",
       "nfev                                                      212\n",
       "njev                                                        0\n",
       "nlu                                                         0\n",
       "status                                                      0\n",
       "message     The solver successfully reached the end of the...\n",
       "success                                                  True\n",
       "dtype: object"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results, details = run_ode_solver(system, slope_func, max_step=3)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXxU9b34/9eZNfu+EZKQgOHDJhBkEURQwKJtL2BptS2t397u2u1aa+t+a6teW6+2Wm2tvb1WfxdrtbZVWi20ilItRVRUQPiwLyEb2fdZz++PMxkymQQGyM77+XjMY2Y+58zMJyHknc9y3m/DNE2EEEKI4cY21B0QQggheiMBSgghxLAkAUoIIcSwJAFKCCHEsOQY6g70N6WUG5gDVAKBIe6OEEKI3tmBMcBWrbWntxNGXYDCCk7/GOpOCCGEiMnFwOu9HRiNAaoSYO3ateTl5Q11X4QQQvSiqqqKNWvWQOh3dm9GY4AKAOTl5VFQUDDUfRFCCHFyfS7FyCYJIYQQw5IEKCGEEMPSkEzxKaXmAn/WWuf0cbwI+DVwIVADfENr/eIgdlEIIcQQG9QRlFLKUEp9EdgAuE5y6tPA+0Am8CXgaaXU+EHoohBCiGFisEdQdwIfAe4CbuvtBKXURGA2cJnW2gu8opR6AfgCcOtgdVQIMTSCwSDl5eW0tbUNdVfEWXA6neTk5JCSknLG7zHYAepRrfUdSqlLTnLOFOCI1rr7T+duYO6A9kwIMSzU1tZiGAZKKWw2WSYfiUzTpKOjg2PHjgGccZAa1H99rXVFDKclAe092tqBhP7vUe86vX7+tuUwO/bXUt/ciZQkEWLwNDY2kpubK8FpBDMMg4SEBMaOHUtNTc0Zv89wvA6qDYjv0ZYAtA5WB6rr2tFHGtBHGgCIdzsoyEmmKDeZwrxkkuKdg9UVIc45gUAAp1P+j40G8fHx+Hy+M379cAxQHwBFSql4rXVHqG1SqH1QVNZFzn13ePzsPdrA3qNWwMpMiaMwL5nCnGTysxNxOuyD1TUhzgmGYQx1F0Q/ONt/x2E3htZaa+A94G6llFspdSmwEnhqsPpQWpjGRdPzGT82lThXdAyva+7k3T3HWff6Af7+5pHB6pYQQpy18vLyoe5CzIZFgFJKrVFKdZ/CWw1MxroG6n+AL2itdwxWfzJT4ylTOXx4QQlfWDGVq5cpFpyfT0FOMnZb5F8EBTnJUa8/XNXM8YYOWbsS4hxQW1vLd7/7XebPn8/MmTP50Ic+xMMPP4zf7x/qrkX58Y9/zBNPPAFARUUFZWVltLS0UF5ejlKK5ubmXl/32c9+lt/85jeD2FPLkEzxaa1fBdK6PV8LrO32/ChwxeD3LJphGGSnx5OdHs+sSTn4/EEqals5Wt3C0aoWCnMjA5Rpmrz2TjnNbV7SktycV5jGeQVpZKbGybSFEKPQt7/9bYqKili/fj0pKSns3buXb3zjG/h8Pq6//vqh7l6E+vp6kpOt31n5+fls27YNgKampqHsVp+GxQhqJHE6bIzLS2HhjLF8avkk0pLdEcePN3TQ3OYFoLHVw1u7qnn6b5q163ezZUcldU0yshJiNHnvvfdYvnx5eCt1aWkpt9xyC/Hx8WzZsoXZs2dHnL9kyRL+/ve/A/Dcc8+xbNky5syZw+rVq9m0aVP4vI0bN7JixQrKyspYtWoVW7duBaxNJI8++ihLly5l3rx5fOtb36K+vh6ALVu2sHz5cu655x5mzZrFpZdeynPPPQfA448/zrp16/jtb3/LV7/61V5HTU8++SSLFi1iyZIl/M///E+vX6/H4+G//uu/WLx4MRdddBF33HEH7e09N173j+G4SWJEs9sNVFE6Byqa8PmD4fbGFg9bd1WzdVc1GSlxTC7OYHJxBnFu+ScQ4lTe3FnFmx9UxXTu1PGZXHpBYUTbxrePsvNA3SlfO3dKHnOnnl6ZniuuuILvfOc7rFixgnnz5jFr1iwWLVrEokWL2LJlS5+vq6+v5/bbb+f555+ntLSUZ599lrvuuov169ezb98+vvnNb/KTn/yEJUuW8MILL3DdddexadMmnn76aV544QUef/xxcnJy+NGPfsT1118fnro7dOgQCxYsYPPmzbz//vt88YtfZPz48fz7v/87WmuSk5O59dZbe12L2rt3Ly+99BLHjh3jc5/7HIWFhSxfvjzinPvuu49du3bx3HPP4Xa7uemmm7jrrru45557Tuv7FgsZQfWzzNR4Lps3ji+smMaHF5RQWpiO0xH5ba5v7uSN9yt45uU9MpoSYoS75557uOmmmzhy5Ajf+973WLBgAddccw179+496etcLhcOh4Pf//73bN++nY997GOsX78ewzB46aWXmD9/PsuWLcNms7Fq1Soee+wx7HY7zzzzDF//+tcpKioiLi6OG2+8ka1bt3Lo0CHA2tr9ve99D7fbzZw5c7jsssv4y1/+EtPXctNNN5GYmMjEiRP5xCc+EfU60zR59tln+e53v0tWVhbJycnccMMN/PGPf8Tr9Z7R9+9k5M/3AeKw2xg/NpXxY1PxB4IcrmxmX3kjhyqa8QWskdWkcRlR61KmacpalRAjiM1m48orr+TKK68kEAiwfft2fvGLX/DFL36Re++9t8/XJSUl8cQTT/DLX/6Sa665hri4OP7f//t/fOUrX6G2tpYxY8ZEnF9WVgZYmxtuvfVW7rjjjvAxh8PBsWPHcDgc5OTkEBcXFz6Wl5cX0849m80W8Zl5eXlRI8D6+no6Ozv5/Oc/H/F7quvzS0pKTvk5p0MC1CBw2G1MKEhjQkEaPn+AvUcb2XmgjiklGVHnvvjPQ8S57Ewdn0luRoIEKyGAuVNPf+qtu0svKIya9usPmzZt4pZbbmHjxo04nU7sdjszZ87kzjvvZPHixTidzogLVU3TDG9IaGpqCq8n+Xw+3njjDb7xjW8we/ZscnNz2b59e8RnPfTQQ6xevZqcnBzuuOMOLr744vCxPXv2UFxczLZt26irq8Pv9+NwWL/eKyoqooJdb4LBIHV1dWRmZoZfl5+fH3FOWloaTqeTZ599lvHjrfzdXq+Xo0ePUlRUdAbfwZOTKb5B5nTYmVKSySeWTiQpITKhe2OLh4MVTew6VM/vX9nL7/6+h92H6wkEZRpQiOFo9uzZOBwObrnlFioqrExuNTU1/OIXv2DWrFkUFRXh9Xp58cUXCQQCPPHEE+EkuPX19XzhC1/gzTffxOl0kpubi2EYpKamcsUVV7B582Zee+01gsEgL7zwAmvXriUtLY1Vq1bxyCOPUFlZSSAQ4LHHHmPNmjV0dnYC0NrayiOPPILP52Pz5s288sorrFixArCmFVtb+07Kc//999Pe3s7OnTt55plnWL16dcRxu93OihUr+O///m8aGhrwer386Ec/4qtf/epAfHtlBDWcHKqM3OpZ29jB3988wps7qyibmMPkkgwcdvmbQojhIiEhgbVr1/Lggw9y1VVX0dLSQnJyMpdeeimPPPIIGRkZ3Hbbbdx3333853/+JytXrmTWrFkAlJSUcOedd3L77bdTU1NDeno6t912G6WlpQA8+OCD3H///Vx//fWUlJTw6KOPkpiYyJe//GX8fj9r1qyhsbGRiRMn8utf/zq8izAxMZGmpiYWLlxIWloaP/7xj5k8eTJgbej4j//4Dz71qU9x3333RXwtdrud/Px8Fi9eHF5buuiii6K+5ltuuYUHHniAlStX0t7ezowZM/jVr36F3d7/GXWM0bZIr5QqBg6+/PLLFBQUDHV3TotpmtQ0dLBjfy17jzbiDwQjjse7HcycmM20CVm4nZJeSYxOu3btCv9CFadny5YtfO1rX+Ott94a6q6E9fXvWV5eztKlSwFKtNaHenutjKCGEcMwyM1IIDejiItm5LNjfx3v7T1Oh8e6Ir3D42fz9kre2V3DxWVjmTQueg1LCCFGCwlQw1Scy8HsybnMKM3ig4P1bNM1tHZYi60eX4CUhJMVJBZCiJFPAtQw53TYmVGazbTxmew50sjbupp4l4MxWYkR55mmSYfHT0KclCkQ4lw1b968YTW9d7YkQI0QdruNySUZqHHpdHr9UdvP9xxp4NV3yimbmMPMidm4ZI1KCDHCyZawEcZmM6JGSf5AkH/tqMLnD/LmB1X8fy/tYvu+WtmeLoQY0SRAjQJtHb6IEVOHx89r28r57frd7CtvlHRKQogRSQLUKJCa5ObqZRNZNqcoohx9Y6uHv24+xO9f2UvF8b4vzhNCiOFI1qBGCZvNYFJxBucVpvH+vlre3l2NxxsAoLq+nT+8uo+SMSnMn55PRkrcKd5NCCGGnoygRhmH3cYslcNnL59MmcqJqAB8sLKZd/ccH8LeCSHOxtGjR4e6C4NKAtQoFed2cNH0fNZcPplJ49IxDAO3y86F08484aYQIlKspdD/93//l9mzZzNv3ryT5sI7mV27dnHVVVed0WsBlFLs2rXrjF8/FGSKb5RLSXSxbO44ZpTm0NrhjdoB2Onx0+H1k54s035CDJSnn36am2++OSr56ulobm6OyIx+LpAR1DkiOz2ekvzUqPbNOyp5eoNm8/bKiArAQojY/exnP+OGG27g2muvpaysjA9/+MPhsu7Lly/nyJEj3Hnnndx0000APPPMMyxfvpw5c+bwhS98IWLqbtu2bVx99dWUlZWxfPlyNmzYQF1dHV/60pdoaWmhrKyM6urqU5Ze/81vfsPFF1/M3LlzefTRRwf3G9JPZAR1Dquub+eDg/WYpsnbu6vZc6SBi2eOpSQ/RepQiWGj/cB7tO99CzMw8KMHw+4koXQ2CeNnnPZrX3rpJX7+85/z0EMP8Ytf/ILvf//7LF26lPXr17NkyRJuueUWli1bxoYNG3jooYf45S9/SWlpKY8//jhf+tKX+POf/0xzczNf+tKXuP7667n66qvZunUrX/nKV/jrX//Kr371q4hEsHfddVefpddfffVVHnnkER5//HFKS0v54Q9/2N/fqkEhI6hzmMNuJaft0tLu5cV/HuTPrx+kqdUzhD0T4oSOg+8NSnACMAM+Og6+d0avPf/887nkkktwOp2sWLGC48ePh2s/dffMM89wzTXXMHXqVFwuF1/+8pdpbW1ly5YtbNy4kdzcXNasWYPD4WD+/Pk89dRTpKZGzn6cqvT6iy++yIoVK5g2bRput5vvfve7Z/Q1DTUZQZ3DMlPjWX3peew6VM/m7ZXhrOmHq5r57YZW5k7NY2ZpNjabjKbE0IkvmTGoI6j4ktMfPQHhSrRAuJptMBg9bV5RUcHPf/5zHnvssXCbz+ejoqKChoaGqOq306ZNi3qPU5Ver62t5bzzzgu3p6SkhOtFjSQSoM5xhmEwpSST8fmp/GtnFTsP1GGaJv5AkH++X8G+o40smV1IVlr8UHdVnKMSxs84oym34SonJ4drrrmGT37yk+G2/fv3k5+fz/r166muro44/ze/+Q3z5s2LaDtV6fWcnJxwhV+AtrY2WlpaBvCrGhgyxScAa1v6JbMK+MSSUrK7BaOahnae+fsequvbT/JqIUSsVq1axeOPP87+/fsxTZN169axcuVKqqurWbx4MdXV1Tz77LMEAgE2b97MQw89RFJSEi6XC6/Xi8fjOWXp9SuvvJJ169axbds2vF4vDzzwwIhMeSYjKBEhJyOBjy+dyDZdw9YPqggETcZkJZKTLiMoIfrDqlWraG5u5rrrrqOmpoaioiIeeeQRiouLAXjsscf4r//6L+69915yc3O5//77KSwsJDMzk8mTJzNv3jx+97vfnbT0+rx587j55pv59re/TXNzM5/4xCdIS0sb2i/8DEjJd9GnhuZONr17jMVlBaQlu4e6O+IcISXfR5ezKfkuU3yiT+kpcaxcNCEqOAWCJi/+8yCHq5qHqGdCiHOBTPGJ0/bunhoOHGviwLEmJo3LYOHMfOJc8qMkhOhfMoISp8UfCLJ9X234+e7D9Ty9QXO4UkZTQoj+JQFKnBaH3cZVyyZSWnhiwbW1w8e61w/wyltH8PgCQ9g7IcRoMqjzMkqpGcCjwHTgAPB5rfXWXs5TwC+AWUAL8KjW+u7B7KvoW0Kck+UXFnNeQSOvvlMevsD3g4P1HKlqYemcIgpzk4e4l0KIkS6mAKWUige+ClwAOIGI1AJa61PmgFdKuYDngZ8Ci4DVwAal1Ditdc/5obXAH4BlwATgdaXUdq31C7H0VwyOCQVpjMlKZNO2Y+wrbwSs0dTzm/YzbXwmC6bnR5SiF0KI0xHrFN9jwF1APNAOtPW4xeISwKm1/qnW2qe1fhrYCVzdy7kqdG8AZujWGePniEGUEOfk8vnFLL9wXMRGiR0H6ni/21qVEEKcrlin+C4HPq21fv4sPmsK0LNa1m7g/F7O/SFwN/ADwA78RGu94Sw+Wwyw0sJ0xmYn8do75ew/1kR6chwzJ2YPdbeEECNYrCMoH7DnLD8rCWv01V07kNDLuSZwQ+g1M4GPKaW+cJafLwZY12jqQ/PGsWxuEQ575I9XMDi6LgoXQgysWEdQPwHuVUp9UWt9/Aw/qw1rirC7BCCi/rFSajZwvdY6P9T0nlLqx8B1wK/P8LPFIDEMg4lF6VHtpmny0uZDpCS6uHDaGJwO2UAqRj6lFHFxcdhs1s9zMBgkIyODq666imuvvXaIe9f/urI/bN26dVCyo8caoK7C2nlXpZRqAbzdD2qtc2J4jw+A63u0TQKe7NFWCLiUUobWuutPbj/WKE6MUPpIAwcrmgA4XNnM0jlFjMlKHOJeCXH2nn766XAqn2AwyD//+U+uvfZapk6dyqJFi4a4dyNbrH/GPgx8Gfg88C3gxh63WGwEDKXU9Uopp1Lqk1hB7489znsDa93pTqWUQylVCnwH+G2MnyOGoSNVJ1L9N7Z6+MOr+3jj/Qr8ASkzL0YPm83GwoULGT9+PHv37gUgEAjw6KOPsnTpUubNm8e3vvUt6uvrw6/ZuHEjK1asoKysjFWrVrF1q3XlTVtbGz/4wQ9YuHAhCxYs4MYbb6S+vp5gMMgll1wSLikPoLWmrKyMtra2k5aC/8Mf/sAnP/lJPvnJTzJv3jx2795NU1MTN998MxdddBGLFy/mgQcewO+3Lh0JBoM88MADzJs3j4suuog//OEPg/WtBGIMUFrrJ7puWAHl+R5tsbyHF7gCa3t5PXArsEprfVwptUYp1Ro6ryZ03qVALbABeBz42Wl+bWIYuWxuEUtmF4a3nZumyTZdw+/+JqU8xKm9dex9Htu6lse2ruWtY+9HHd989O3w8fereu7Fgk2HtoSP7zq+N+r4y/tf7/O9T4ff7+fFF19k//79zJ07F4Ann3ySF154gccff5zXXnuNjIwMrr/emkzau3cv3/zmN/nmN7/J22+/zec+9zmuu+46Ojo6uOOOO9i7dy9/+tOf2LBhAx6PhxtvvBGbzcaKFStYt25d+HPXrVvHZZddRmJiIvfddx87duzgueee469//St1dXXcdddd4XO3bdvGddddx8svv4xSiu9973u0tbWxfv16nn32Wd58801++ctfAvDb3/6WP//5zzz77LOsX7+enTt3ntX353TFfKGuUupa4BYgP/S8BnhQa31vrO+htd4BLOylfS3WtU9dz/8FXBzr+4rhr6swYmFuMq+8dZSj1daIqqGlk9+/spdZKoe5U3Kx22VtSowsn/70p7Hb7Xi9Xvx+P/Pnz+fnP/85559vbVB+5pln+MY3vkFRUREAN954I7Nnz+bQoUO89NJLzJ8/n2XLlgFWKY5x48bh9/tZv349a9euJSsrC4Dbb7+dhQsXUl1dzapVq7jyyitpbW0lMTGRv/zlL9x9993hUvBPPvlk+HU33HADH/nIR/j+978PQHp6enjqsba2lo0bN/LGG2+QlJREUlISX/va17jlllv42te+xosvvshnPvOZcN//4z/+g1dffXWwvrUxX6j7HeB2rK3fr2Ndn3QRcJNSqkNr/eDAdVGMJskJLlZcPJ6dB+p44/0KfP4gpmny9u5qDlU0sXROETkZvW3sFGJ4euqpp5g8eTI1NTXccMMNuN1uFixYED5eUVHBrbfeyh133BFu616avWeJ97KyMmpqavD5fOTn54fbs7OzcblcVFZWMnPmTJRS/P3vf2fs2LEEg0EuvPDCU5aCB6uib/e+AVx++eXhNtM08fl8eDweamtrycvLCx8b7BJGsY6gvgZ8VWvdfR3oDaXUYawLeCVAiZgZhsG0CVkU5iaz8e2jlNdYGznrmjt5f18ty+YWDXEPxXAze+x0Zo+d3ufx+YUXML/wgj6PLyqex6LieX0eXzphIUsnRE3unJacnBx+9rOfsXLlSu65555wQMrJyeGOO+7g4otPTArt2bOH4uJi3n33XbZv3x7xPg899BCrV6/G5XJx7NgxsrOt6wmrq6vxer1kZmYC1mjrpZdeIi8vjxUrVmCz2U5ZCn7btm0RgSsnJwebzcY//vEP4uOtTdatra3U1dXhdrvJyckJB7auPgymWOdTsoGonHnA24BUBRRnJDXJzcpFE1hUNhan3Ua828HCGfmnfqEQw1RaWhp33303Tz31FK+//jpgBZJHHnmEyspKAoEAjz32GGvWrKGzs5MrrriCzZs389prrxEMBnnhhRdYu3YtaWlprFixgvvvv5+6ujpaW1u5++67KSsro7CwEICPfOQjvP3222zYsIFVq1YBnLIUfE95eXnMnTuXe++9l7a2NlpbW7n55pu57bbbAPjYxz7G//3f/7F//37a2tr46U9/OgjfxRNiDVA7gE/00n41VjYIIc6IYRhMPy+bqy9TVrokd+Sg3usL4PNLhnQxcixcuJCPfexj3H777bS2tvLlL3+Z+fPns2bNGubMmcMrr7zCr3/9a1JSUhg/fjwPPvgg999/P7Nnz+aJJ57g0UcfJTExkZtvvpkJEyawcuVKLrnkEux2Ow8//HD4c1JTU1mwYAH5+flMmDAh3H7LLbeQl5fHypUrWbBgAYcOHQqXgu/N/fffT2trK5dddhlLlizBMIxwILryyiv59Kc/zWc/+1mWLFnCpEmTBvab10NMJd+VUh8C/gK8BGwONc/HSoH0Ma31nwesh6dJSr6PLhvftjZUXHpBoWRIP0dIyffRZcBLvofy4C0FPMBngY8DzcCc4RScxOhSXtPCzgN1NLd5eX7Tfl556widXv9Qd0sIMUhi3mautd4EbBrAvggRwesL4nbZ8XitKb6uelOLZxVQkp86xL0TQgy0PgOUUuoZ4Ita6+bQ4z7FUg9KiNM1fmwqeZmTeG3bMfZ3qzf1lzcOUlqYxsUzx5IQ5xziXgohBsrJRlBtWFnFux4LMegS4pxcMb+YfeWNbNp2jPZOKyXj3qONHK5q4cJpeUwbn4XNZpzinYQQI02fAUpr/e/dnv4nUK61jkicppSyY5XDEGJAnVeQRkF2Eq+/d4zdhxsAa4ffpm3H2H2ogRUXj4/aAShGLtM0I67XESNTLJvwTibWbeYHgcxe2ouAf5xVD4SIUZzbwbK541i5aAJpSe5wu9tlx+2S0vKjhd1ux+eT4gWjQUdHB07nmU/Dn2wN6vNA19VdBrBeKdVzC1UecOiMP12IM1CYm8wnP6TYpmt4d89xFpcVyF/bo0haWhrV1dWMHTs2XGdJjCymadLR0cGxY8fIzc094/c52ZzI77CyRBjAbOBlIosLmqHnvz/jTxfiDDnsNuZMyWNGaXY4Q3qXQCDIq++UM0vlkJ4SN0Q9FGcqKyuL8vJytNZD3RVxFpxOJ7m5uWdV2PBka1BtwA8AlFKHgKe11p4z/iQhBkDP4ASwbc9xdh2qRx9poGxiDrMn5+B0yBTgSGGz2cLZs8W5LaZVZa31E0qpMqXUFKxigmCNrNzABVrrrwxUB4U4HR5fgLd3Wwktg0ErS/quQ/VcOC2PSeMyZLefECNIrOU2bgV+iDWllwg0AV1XSr44MF0T4vS5nXZWX1rKq++UU1VnXR3R3unjlbeO8v6+Wi6ani8pk4QYIWJdgfwKcKPWOgWoxCrVPhb4F71nORdiyGSlxbP60vNYOruIpPgTO4hqGzt4ftN+1v3jAPXNnUPYQyFELGINUHnAc6HH7wLztdZVwHexcvMJMawYhsHkkgzWXD6JuVPzcHar1Hu4qpmnN2g2b68Ywh4KIU4l1gB1nBPXQe0BZoQeHyNUAl6I4cjpsDN3Sh5rrpjMlJKM8Hb0oGkSLxf2CjGsxRqgngceU0rNBDYC1yilFgPfBg4PVOeE6C9J8U6WzC7i6mUTKchJJjXJzfkTsqLOCwbP7sp3IUT/ifVPyO8ADwDTgLXAaqzrotqATw1M14Tof1lp8axcNJ4Ojx+7PfLvs8raNv725mEumJTLpHHpUceFEIMr1m3m7ZzIKgHwOaXU14FOrbUU6BEjimEYvWZBf/ODKprbvGx8+yhv7apm9mQJVEIMpZOlOvpwLG+glEJrLVvNxYjW3umjtrEj/LylXQKVEEPtZCOoWCvlmpy4eFeIESkhzsk1H57M9v11bNM1dHisiYHugeqCSTlMLs6QQCXEIDlZqiP5XyjOKU6HnVkqh/MnZPYaqF59p5y3dlVTpnKYUZo9xL0VYvSLNZNEwsmOh9aohBgVThaoWjt81NTLj7sQgyHWXXytnKiu2xuZ4hOjTvdAtWN/He+EAtX0XkZPDS2dpCW5peyHEP0o1gB1aS+vmwDcAHyvX3skxDDjdNgpUzmcf14Whyubyc2InFDw+gI8+/JekuKdTD8vCzUuA6dDZsiFOFuxbjN/rZfml5VS+4D7gD/1a6+EGIYcdhsTCtKi2ncdqsfrC1DvC/DqO+X8a0cVU0oymFKSSVqyu5d3EkLE4mxzvVQAU/qjI0KMVIGAidNhw+cPAtDp9fOOruEdXcOYzEQmFWdwXmEa7l5qVwkh+hbrJonerolKxZrie69feyTECDNrUg5TJ2Sy+2A97+07TnObN3yssq6Nyro2/vHuMUryUymbmE1Oxkn3HAkhQmIdQfV2TZQXq9TG1/uvO0KMTG6nnRkTs611qqpmPjhYz+HKZoKmtbfIHwiy92gD5xWkSoASIkaxrkH1y4qvUmoG8ChWPakDwLO5tGsAACAASURBVOe11lH1pJRSycDPgBVYuwd/D3xda+3rj34IMVBsNoOS/FRK8lNp7/Sx90gjuw/Xc7yxgziXg+IxKRHn+/wBdh6ooyQ/ldQkWa8SoruY16CUUnZgGVbC2ABWXajXtNYxpX9WSrmwsqL/FFiElXB2g1JqnNa6ucfp/ws4gWIgDngJuBG4J9b+CjHUEuKczJiYzYyJ2dQ2dtDU6onKQnG4qoXX36vg9fcqyEqLZ/zYVMbnp5KZGidb1sU5L9Y1qGLgr1gB4zBWmY4iYLtS6iNa6+oY3uYSwKm1/mno+dOhhLNXA7/q9lljgJXA2FDgalZKrUSutRIjWFZaPFlp8VHtB481hR/XNnZQ29jBmzurSEl0WcFqbCp5GYnYbBKsxLkn1hHUL4GDwCKtdQ2AUioPeBJ4BPh4DO8xBdjVo203cH6PtjLgCLBGKfVNrJHU/wG3x9hXIUaMcWNS8PgCHK1uIdCtFlVzm5d39xzn3T3HSYhzUpKfwtSSTFm/EueUWAPURcDcruAEoLWuUkrdALwR43skAT1zxLQDPf/HZWCN1KZhrVXlAOuAFmSKT4wyE4vSmViUjtcX4HBVMweONXO4qhmvLxA+p73Tx84DdRTkJEUFqGDQlNGVGLViDVAHsUZAH/RoL8C6FioWbUDPOY4ErDRK3XmwpvNu0Fq3Aq1KqQeAa5EAJUYpl9NOaWE6pYXpBAJByo+3cuBYEwcrmmnvtPYGjc1OinhNIGjymz/vJDM1noKcJApzk8lOi5eAJUaNWAPUT4GfK6UmAq8DfmAWcBtWKfjwdVInqQ31AXB9j7ZJWNOE3e0O3acBXRP0Z3tBsRAjht1uY1xeCuPyUlhcZlJd305NfXtUkcWa+nY6PH7Ka1oor2nhXzsqcbvs5GUkkpuZQF5GAjkZCcS55L+PGJli/cnt2sRwVy/Hbuv2+GS1oTYChlLqeuBhrF1804E/dj9Ja71dKfUW8BOl1GeBLKzA9iuEOMfYbAZjshIZk5UYdaymITqrusdrTRUerrI2xhqGQXqym3FjUrhoev6A91eI/jRo10Fprb1KqSuwroP6AXAIWKW1Pq6UWgP8UmvdNYfxYeBBrGulbFjbzn8a/a5CnLtmlGYzoSDNGkFVt3C0upW2zshLBU3TpL65k9REV9TrK4630tzuJSs1nvRktxRiFMPOaY39lVJLgalYQWMX8LLW2h/r67XWO4CFvbSvBdZ2e34c+PTp9E2Ic1FSvJNJ4zKYNC4D0zRpavVSVd9GdV07VfVt1DV2EjRNcjOjR2AfHKxn9+F6AGyhkVZmWjxZqfFkpsaRmRZPYpxDrscSQybW66DysDKWz8Ia+RjAOGC3UmpZ9919QoihYRgGaclu0pLdTBqXAYDPH+R4QztJCdEjqLqmjvDjoGlS19xJXXMne2gIt8e5HKQnu5k7NY/C3OSB/yKE6CbWEdSDWBsjSrTWxwCUUvnAU8ADwGcGpntCiLPhdNjI77H7r8uEgjRSEl3UNnXS1Orp9ZxOr5/KOn84p2B3f9i4F5vNRlqSi7RkNymJbpITXCQnOHG77DLyEmct1gB1OXBpV3AC0FpXhK6D+tuA9EwIMaBmT84NP/b6AtQ3d1Lb2EFdUyd1TR3UNnWGr8dK65EnMBAIUlnXjmmalPcyf+J02EhOcJGU4CQlwcWF08YQ55bdhOL0xPoT00nvJd9PtmtPCDFCuJx28jITyeu2VmWaJq0dPhpbPCT3mCJsavNi9jKq6uLzB6lv7qS+uROAi2ZE7iBsbPHw7Ct7SIxzkhDnJCHOEXrsIDHeuo93W7c4l0Ou7TpHxRqgNgAPKKU+2ZV3TymVC9wPrB+ozgkhho5hGKEpu+j1q7QkN2sun0Rji4emVg+NLR5a2n20tHtpafeGizeCtY7ldET+HdvW6cPjDeDxBsJBrC9J8U4+99GpEW3HGzrYc7SBeJcDt8tu3ZyhW+ixy2mXwDbCxRqgbgReAQ4rpQ6H2sYB7yO77YQ459hsBunJcaQnx0UdM00TjzdAS7uP1o7IYNWlrSP2yjnxvUwN1jZ2sC2GvVkup51xecksv7A4ov1wZTNHqltwOWw4nXbr3mHD5bTjsFuPnQ4bDrsNt8suFzsPkVivg6pSSk3HWouaAnQAu7TWfx/IzgkhRh7DMIhzO4hzO8hOj87gDlBamEZBThLtnX7aOn10hO7bO0L3nT7aPX46PYFeA1SHN7arW7y+AMFg9FRkRW0b7+09HtN7TBufySUXFEa0vfF+BQePNeFw2LDbDBx2G3a7gd1mBTWH3cAeui/JT41KU3W4spkOrx+7zXqN3WZgsxnhe5tx4nlivBOXM3IE6g8EMQwDm8Go3owS858FWmu/Umov4MaqB7V/wHolhBjVDMMIrT05yYpK0Rmpt7WusdlJzD9/DB0ePx5vAK8vgMdnTRl6Qo+9viCmaeJ2RS+T+/yBqLa+9JyeBGht99HYx87HnhLinFEB6h1dw7HjPdOQ9u5D88YxsSg9ou3pDTr8+TbDwDCsUW04uHVr+9C8cRFriwB/fHUfgaAZDnBG6HzD6Hq/rucGi2aOJTH+RJotnz/IP94tZ+GMsVGBs7/Feh1UBvAsVk0nL9Z1UE6l1Dpgjda6bcB6KIQ4p/U2QsjNSCD3FKVHTNPE6w/2GuC6tth7fUG8/gA+fxCvL4gv9NgfCOL3B/EFgsTHRf+a9Aeipy374rBH9z/Qy6iuL7Zevv7u2/6Dpglm3+/Z2yUC1fXtMX8NC84fE/E8EAzywcF65p8/8KmzTqceVCowS2v9HoBSahZWfrwHgS8OTPeEEOLMGIaBu4+/8MdmJ0WNak7HpRcUsGD6GAIBE38gSCBohoNa1+NAwCQQDEaNXgDG5SWTmugiEDRDtyDBoEkw9Dz8uI8RYNcoqbfg09u5PfU27dmXnn8gdH3kYMwsxhqglgOLu4ITgNb6HaXUV7F2+EmAEkKcM7qmJ8/UnCl5Z/X5n7liMmCNEk3TGiV1BbWgaRI0rSBkmmbE9FyX1UtKw8e7Rl5dI7Fg6D273jveHRkgnQ4bl15QiHMQcjfGGqCOA5m9tLuwCgkKIYQYZF1rRTaM07oi9VTToyfjsNuYOr63cND/Yg1Qt2LVfbqNyHpQ92PViZrSdaLWumdRQyGEEOK0xRqgngrd/x8nMkp0zUDeA9wdei6ZJYQQQvSLWANUyYD2QgghhOgh1gt1D5/6LCGEEKL/SAlNIYQQw5IEKCGEEMNSnwFKKTVNKSUBTAghxJA42RrUP7ESw5YrpV4BPqa1bhycbgkhxLnLuoA2QMAMYrfZcdgiN0fXdzTS4eskEAyQlZhBgjMyn+GB+iO0eFsJmkHOyygm2R2ZNWNb5Q6aO63jF4ydTkqP468d/Bct3jaCZpBLS+ZHvX7d7r+zvHQxLvuZX6wci5MFKA/weaXUa1g5+BYrpRp6O1FrvWkA+iaEEMNSIBjAG/DhC/px2hzEOyPLjlS21FDTVosv4GdsSh5jknMijr9ftYuDDUfxB/2UjZnG+IyiiOOvHHiD/fXW3rQl4y/ivMziiONby9/jcGM5AB86bzHF6QURxz84vpeK5ioAshMzowLMoYZyjrfVATAlpzQqQFW31dLY0QSALxidOb4r+A20kwWoO4H7gO9jXd/0xz7Ok2ufhBAjSiAYwOP30hnw4LA5on5BlzdXsrf2IN6Aj6K0fCZnl0Yc31a5k3cqtgNQlj+NOWNnRBw/0nSM9yqtnAV2my0qQLV4W6lutcp9tPvao/pn7zZi8gejM6/bjBOrL0Gzt+MnEuUFg9GBJPL10Xn5uh/vLdmugXHSisr9pc8ApbV+WCn1CJCAlc5oAnDqCmFCCDHEGjubOdhwhA6fh7S4ZKbkTIw4vrfuEJsO/QuAiVnjuaRkfsTx5s4W9tYdBLCmz7Ij37/71Ja/lxGG03biV2uglwBjN7oHoOgA4rDZsdvs2A1br0lZMxLS8AS82AyDOGd00ciS9ELS41OxGbao4AswI28yHf4SbIaN1LjkqOMLx80hEAxgGEavr//wxEtxO6IrLfe3k14HpbU2gTalVAlwBIgDSrE2V+zXWksePiHEgDFNE2/AF/XLsLa9nn8dfYdWbzvpcaksL10ccbyps5mt5VZu64LUMVEBqvv7dfqj6zq57CeOewPeqONxDjdxDjcOuwOnLXodJicpi+l5k3HYHOT3GD0BTMkupTitAIfNQaIrOi/ewnFzWThublR7lwvyz+/zGBA14uupOL3wpMfzkrJPejw1LuWkx/tLrJkkyoF7gW9hJYgF8CqlfgN8XWsdW3lLIYToxjRNPAEvcQ53RHtTZzN/2fMK7d4OUuOS+cS0j/Z4IVQ0VwNgN6I3G8d32zTQ4euMOh7ncBPnjCPO4SLRGR0gcpKyWFxyIS67i2R3dLmMiVnjmZg1vs+vqyBlDAUpY/o8nhKXTEovIxcRKdYAdTfwGeAarGSxBnAR8N/AHaGbEEJEMU0zqqZQp9/DC7v/RounFafNwTVlH4847rK7aPVYdVDbvNFrNN1HHR29jICSXYlMz5tMvDOOFHd0IBiTnMM1M1f32ecUd1KvU1ticMUaoK4Bvqi1frFb2zNKqRbgMSRACSG6MU2TTYe20NDZRFNnC5+ZcWXEwr/b7qLZ00IwGAzviOu+rhPncGOz2QgGgxiGDX8wELHVOs7h5oqJl5LoSojaYg0Q74zjwsJZA/tFigEXa4BKAvb10n4AyOq/7gghRorjbXWhWz3zCssipukMw+BYS1V4FNTkaSEjPi3ieIo7icaOZlwOF+2+jogAZRgGV037NxIccTjs0b+mDMOgMHXgS46LoRVrgNoKfA1rDaq7rwNv92uPhBDDij9gLTH3DBT/OPwmtW31AEzIHBe15pIelxoOUI0dzREBCmD5eYuJc8T1uRtMpthErAHqe8CrSqlLgH+F2i4EioHL+79bQoih9n7VLnTtfho6m7mk+MKoTQHZCZnhAFXb1hAVoMrGTGV63mTS41KjLmSFwdsJJkauWMttvKWUKgO+DEwFOoB1wCNa68oB7J8QYgCZpklDZxMGBunxqRHHOv0eGkLZBGrbG5jY47VjU/LwBrxkJ2b2umMtr5ft1UKcjlhHUGit9wI3DmBfhBCD6GDDUTYd2oLH7+G8zBKWjF8QcXxMcg7vVu4Ew8DTy0658RlFUSl6hOhPMQeo/qCUmgE8CkzH2mDxea311pOc78SaUlyntf7+oHRSiFHG4/fS5GkhJzEzoj3RGR8OPFWt0Uli8pKy+eikZWQnZOAc4KSgQvRm0MppKKVcwPPA74A0rGurNiilTjYRfRcwcxC6J8So0+Hr5IXdG3jy3d/z172vRuVOy0rIwGFzEOeMIyshIyrnm9PuJD85V4KTGDKDOYK6BHBqrX8aev60UurrwNXAr3qeHNqQcRmwfrA6KMRIFQgGsBm2iAti4xxuGjtbME2TTl8nte31ZHcbRdlsNq46/6MkOhOiLqQVYjgYzIKEU4BdPdp2A1FJpZRS6VhB6xogOhGWEAKA8qZKNh74J0+++xz1HZHl2gzDsDYvGAbZiZm9lk1IciVKcBLDVkwjKKVUNnAbcAHgxEp1FKa17jur4QlJQM+cJe1Y2dJ7ehT4udZ6h1Iqli4KcU7StfvDdYMONBwhMyE94vicsdOZXzir123eQgx3sU7x/RqYB/wf0HyGn9UG9MxJkgC0dm9QSn0OKzvFTxFCEDSDVDRXY2JGZU+YkDEuHKC6CtB117NQnRAjSawBahGwUmv92ll81gfA9T3aJgFP9mj7FDAXaAiNnhKBy5VSs7XWPVIaCzG6VbUe52/7NtHh6yQrMSMqQBWk5jNzzFTGpxdFjZ6EGOliDVANQNNZftZGwFBKXQ88DKzG2m4eUalXa728+3Ol1J+Ad2WbuTgXpbmT8YTqEdW21dPU2RyRgcFhszO3QDa6itEp1gB1B/BwKLjspcfGBa11dD78HrTWXqXUFVjrSz8ADgGrtNbHlVJrgF9qrWU+QpyTGjqa0LX7KRszLSI3XZwzjuK0QipbqpmQURxRiVWI0S7WAPUA1rVL/+rjeEz/a7TWO4CFvbSvBdb28ZpVMfZRiBFp06Et7D5uFQtIciUyLTdyY9BFRbNx213YbIO56VaIoRdrgPr4qU8RQpyJ7mtHu2v3MzVnYsTWb9mBJ85VsSaLfQ1AKRUPlGJdP7Vfa90ygH0TYlRp6Giitr2e0sySiPbSjGK2lr9Lfkoek7PPG6LeCTH8xHodlB24B6seVNd1UF6l1G+Ar2uto68AFEIA4A34eOXAGxxpPIbNZqMgZUzEqMjlcLFmxpWSUkiIHmKd1L4b+AxWZofC0O0a4MNIuXchTsppc9Du6wAgGAyys2ZP9DkSnISIEusa1DXAF7XWL3Zre0Yp1QI8hgQpIQDwBwP4Az7iuo2QDMNgRt4UXt7/OuPSCihIja6dJISIFmuASgL29dJ+ACvrgxDntEAwgK49wLbKHeQlZbN0QuRm1ZL0Qq46/99IkyqyQsQs1im+rcDXemn/OvB2/3VHiJGprqOR1w+/SZu3nf0NR6ISt9oMmwQnIU5TrCOo7wGvhkpgdF0LdSFQDFze/90SYmTJScykKG0sRxqPEe9w0+ppIyM+bai7JcSIFus287eUUrOAL2OVzegA1gGPaK0rB7B/QgwrpmlysOEIia4EcpOyI47NHjudMck5TMkulU0PQvSDmAsWaq33AN8ZwL4IMazVttXz6qF/Ud/eQG5SNismXRZxQW1WQgZZCRlD2EMhRpc+A5RS6k1guda6QSm1FTD7OjfGelBCjGhuh4vGDitncnXrccqbK6Oyiwsh+s/JRlB/ATyhx38ehL4IMawlu5OYlD2BPbUHmZo7MaJ8uhCi//UZoLTWd3Z7uhHYrLX2dT9HKeXGulhXiFHDHwywvXoXGfFpjEsriDh2Qf50ZuWfT4KzZ+1NIUR/i3UNaiOQBxzv0T4eeIroSrlCjEjH2+r42/5/0OppI9mdREHKGOy2E8n6JXGrEIPnZGtQ1wJdoygD+EAp1XMdKgnYNkB9E2LQJbuT8AasiYIWTyt76g4wObt0iHslxLnpZCOoXwFtWBfz/i/wQyKr6ppAK/DygPVOiEEW53AzO/983qnYweyxM1BZE4a6S0Kcs062BuUHngRQSh0E/gkka63rQ21zgG2SyVyMREEzyAc1e3DanVFBaEr2REozx0dUthVCDL5YUx0dBzRwc7e2dcB2pZT8iSlGlBZPK3/84K/888jbbD7ydjjTeBebzSbBSYhhINYA9TPgH5xYkwIoAbYAD/d3p4QYSPHOeHxBa+DvDfh4v2rXEPdICNGbWAPUXOAHWuvWrgatdQdwF3DRQHRMiIHisNlZWDQHu83OnIIZzBk7Y6i7JIToRazbzOuBacD+Hu0TASn7Loat+o5GalprmdSjlHpB6hg+PX2VbBsXYhiLNUD9D/CYUqoQeAtrB98s4Hbg1wPUNyHOWCAY4K2K98PTdzlJWVHZxSU4CTG8xRqg7g6dezvQlcK5BvgJcN8A9EuIs2IYBpUtNZimdeneP4+8zUfV0iHulRDidMRabiMI/Cfwn0qpLMCrtW4e0J4JcRZsho3FxRfy3M4XyUnKYkHRBUPdJSHEaYq53IZSqgyrFpQ99NwA3MAFWuuvDEz3hDg10zQ53l5PTo/krenxqayaspzM+PSIshhCiJEhpgCllLoVK5NEK5CIlVEiNXT4xYHpmhCn1uZt5/XDWzncWM6H1RIKUsZEHJf6TEKMXLFuM/8KcKPWOgWoBKYDY7HKv28doL4JcUpbj73H4cZyADYd2oIv4DvFK4QQI0WsASoPeC70+F1gvta6Cvgu8NmB6JgQsZhbMBNXKOtDQcqYvqtqCiFGnFjXoI4DmcAhYA8wA3gWOAZISVExKIJmEAMjYj0pwRnPJcUX4nK4yE/OHcLeCSH6W6wjqOexroOaiVUb6hql1GLg28DhgeqcEF3qOxr506717Kk7EHWsOL1QgpMQo1CsI6jvAA9gZZNYC6zGKrPRCnx6YLomhKW8uZK/7nmVoBlk85G3KUgZQ6IrYai7JYQYYLEGqNXArVrrutDzzymlvg50SrkNMdDyErNJcifQ3NmK3wxQ3VrL+Iyioe6WEGKAxRqgHgLeBLoCFN0Tx8ZKKTUDeBRrF+AB4PNa66hdgEqpC7CyVEwHmrFSLf1Qay1r4Ocgh93BouILebP8XRYXX0h6fOqpXySEGPFiXYPaAlx5Nh+klHJhrWX9DkjDSp+0QSmV0uO8BOAvwDNYGzOWAp8DvnQ2ny9GhprWWnYd3xvVnp+cy8pJH5LgJMQ5JNYRVBC4Ryl1G3AQiKjwprWeG8N7XAI4tdY/DT1/OjRNeDVWefkuhcBmrXVXnam9Sqk/AQuBx2LsrxhhAsEAW4+9x/vVu7FhkJuYTUZCZHJXyQYhxLkl1gC1JXQ7G1OAnpXhdgPnd2/QWmu6jdZCI68rkOA0qhmGQXVrLZgmQUzeOPIW/zZp2VB3SwgxhPoMUEopWyhJLFrrO/s67zQkAe092tqBPrdjKaXcwFOh8x7thz6IYcpm2FhUPI/nPniRMUk5LCqeN9RdEkIMsZOtQfmUUjndG5RSi0JB40y0AfE92hKwtqpHUUrlAa8AOcCyUAVfMUrUttdHtaXHp3Ll5Mv58MQlJLuThqBXQojh5GQBqrcJ/z9j5eA7Ex8AqkfbpFB7BKXUFKwcf/uwglPDGX6mGGY6fJ28vP91/rDzJcqbKqOOZyZI5nEhhCXmchshZ/ObYyNgKKWuBx7GurZqOvDH7icppdKBDcDTWuvvnMXniWFo67H32F9vJR/ZdHgLn5j6EZx25xD3SggxHMW6zfysaa29WJsdVgP1wK3AKq31caXUGqVU11TfZ7FGadcqpVq73X47WH0VA2f22Om4HdYs8ZikHIKmXNomhOjd6Y6gzorWegfWdvGe7WuxUiihtX4I68JgMcJ1lVvvmdx1UfE87IaNorQznS0WQpwLThWgPtdtZNN1/meUUrXdT9Ja/7zfeyZGtLr2Bv5x+E0mZo5nSk5pxLGS9MIh6pUQYiQ5WYA6Alzbo60K+PcebSYgAUqElTdV8tLejZimSWNHE8XpBSQ4e27gFEKIk+szQGmtiwexH2IUyUvOIdmdRHNnC/5ggKqW45LcVQhx2gZ1DUqcGxw2OwvHzWFbxU4uHjeHNMmfJ4Q4AxKgxBkLmkE+qNmDN+BjVn5ExioKUsYwNjlPrmkSQpwxCVDijHj8Xv6y52Vq2+oxDIPitEJJ7iqE6FeDdh2UGF1cdicuuwuwtpO/VxWVEEQIIc6KBChxRgzDYEHhBTjsDmaPnSHJXYUQ/U6m+MRJmabJ4cZyypurWDhuTsSxjIQ0Pj19FXGOM80fLIQQfZMAJfoUNIO8tGcjx5qrABiXNpbC1PyIcyQ4CSEGikzxiT7ZDBtJrsTw83cqdgxhb4QQ5xoZQYmTmlMwg0ON5ZyXOY4L8qcPdXeEEOcQCVCCQDCArj3AvvpDfHTiUmy2EwPrBGc8nzp/BS6Hawh7KIQ4F0mAOseZpsnzuzdQ22ZVuN1duz8quasEJyHEUJA1qHOcYRiMTx8Xfr67dl+4TIYQQgwlGUGdQ0zTpMXbRoo7KaJ9ak4pe+sOMDFrPFOzJ0oGCCHEsCAB6hxgmiblzZW8XbGdxo4mPjV9Fe5u03ZOu5OPT/2IBCYhxLAiU3zniDfL36WmtRZvwNdrWiIJTkKI4UYC1DnAMAxm5E0BiNihJ4QQw5lM8Y0iLZ5WdtRoPH4vl5TMjzg2Pr2I+jGNTMkuJcmd2Mc7CCHE8CEBapRo87bz9PYXrB14hsGsMdNIiUsOH7fZbMwtmDmEPRRCiNMj8z2jRKIrgbEpedYT02Rv/aEh7Y8QQpwtGUGNIKZpUtNWi649QH5yLudlFkccn543maAZZHru5KikrkIIMdJIgBpBdh3fx+uH3wSgsbM5KkAVpIyhIGXMEPRMCCH6nwSoYco0zait38XpBbxxZCumaVLVUkNzZ0vEOtNwY5pBzIAf0++DYAAzGIBgENMM3QeDYAYjjwHhr9oIzUCfaIDQ98Sw2cDmwLDZMOwOMOwYdjuGzQ52h3Vvs8v2eSFGMAlQw8zRpgoONhylvLmSj0/9CC67M3wswRnPhIxi7IYNlT2B5B4ZIfqTaZqYAR+m14Pp9xD0dmL6vAR9nZg+j3Xv9Vjn+H0n7v3ebm3+AetfrAybHcPhwrA7MJxuDIcTw+602hzO0M2FzeGyjjvd2EL34cd2+W8ixFCQ/3nDzDsVO6huPQ7A/vpDTM6OTNy6ZPyCM3pf0zStwOLtwPR2EPR0EvR1EPR0YHo7CXo7rGOe0GNfJ4yCnHxmMIDp7bCedLSc0XsYNnsoWMVhuEJByxUXeh6HzRWH4bTuba74cHCT0ZsQZ0cC1BDw+L0caTpGoiuB/OTciGMl6YXhAHW48VhUgOrONM0TwcXTTtDTQdDbjukJPfd2hNqs50MRcKzRigPD5gCbDcOwW/e2bveGzZqyM2zWFF73foYfm92aTTCDmIFAeHrQDPhDj/3W1GHAb00b9gMzGMD0tFvfw9i/citwhQKYzRmH4Y7HFgpk4cDmisPmjMfmirOmJiWoCREmAWqQHWoo52/7N2GaJsXphVEBqjg1n9a2BsYlZJJlj6OzXJ8INJ720OjHCkRBTycwcEGna1rM5ow7Md3lOvHcmjJzYXM4we6MvO+aShvCX7imaULQH5p6tKYfg13TkH4fZiB07/Nax3wezNDNmsr0Ws/N4Jl8OkFvB3g7iDVMGjZ7KHDFR4zKrCDmxnDFR/4buOKG/HssxECSADUATNOkobOJho4mSpJyT0yt+TwktDXga67FDAQ40FBFbUMzDp/vxBSb30fXmKm5n/tlOJzYXAmhX3rx2NxdBr3h/AAADyNJREFUvwzjQ8+7HseNirUXwzDAbgVK3Gf2HqZpQsAfDl5da29BX9fUaGhNrutxaHrU9PtO/7OCAczONoKdbbG/yDBOTDU63VZQc0avp4XX1RxubK7QY5v9tPsoxGAa2b+BBoFpmtZf3D3/yvZ7QpsGPKHjHoI+D15PG3+q13QEPNiCJqvdY3AYtu5viMtbRxw2CuzxeGqOYBpnfr204XRbU0juBCvYuBOsQONOCAWdrvb4ER9whoJhGOBwYnc4IT72TSnW2ldnKHCFgpmvK4h1rfl1httNb+eZTUmapz9S6xLeQOJwYThd3TaLuCLb7c7wSLlrU4nVZt2wydSkGBiD+htLKTUDeBSYDhyA/7+9u4+Sq67vOP6+M/uUzdNCnprShBA035iKWBBEQaFSjg/0HESwKBHjAVq1oJZWTosgD9oIVhGKqFHIKQaI8RTCk1pLqymntmBj5DnkWyVQU5otkBA4ybLZp+kf33tn7052kk3YnZ3ufF7nLDP33t+985tfmPnO73d/D5zr7uuHSTcfWAkcBzwPfMrdf1SLPPa+1Mku/zn9Xa+kTUF9VDajvVrqZ/tAD9tKPbyuOIX2ZOgv0ea+XrpK/fQDLw708FvFtvKxJEl4X8scmqsFpSQZvNme1m6ygBM1nHaS1raoCbVO0q/gOpUUiiRtkym0jWzew3JNrad7aO0sC3K9WTDbPRjUene/pp6S5Q4kWSeSA5bEfcamFpJCU9pTsin9a457a+V96WM2HKAY9yeTYhMUi0O3C9mwgeLgvUoFwoZSswBlZi3APcD1wDuBM4D7zexQd69szVoDPAicCpwA3G1mb3b3zWOdz65fb6B3+9by9nDjkR7s3U7nQDcAHUkz84vtQ47PKbTS3d/P7EIbzU0tFNumDmmCmZQFoNbc/Yas6U29vxpSvqZWZORj20pZ82MWvHqHDgcYrOXvLtfys+Oj12mmVL7PN7aS6ExTTMe/ZR1vClkwKww+JhXbhSJJUhjaUSfdLnfSyR2DJP5NsuP5tEk6Hi+7Bkm5g08c2/N59jdkX/YaUlUta1AnAc3ufn26vcbMLgTOAm7KEpnZIuAtwCnu3gP81MzuBc4DLh3rTLbMms8znb/i6f6d7Cj1sqg4hTe0HkzS3EKhqZWkuYU5PS282L0NkgI7O+YxedaiCD5NLRRa2jixUKS1bUoEHNVwZAwlxSaKxSYYYU0tU+5A0tvDQF/PkGbsUl8vA2kTdrlDSX58W+X2KPWWHEGu0wHd/WPYNWgcJEka6EiDF5SDVxrIykEOcmmzY5TTJBVpy48k6dPKY9nz8n+GBs18morrNU2bwaSFR1LILX462moZoJYAT1Xs2wQcMUy637j7rop0x45h3somLTiCYhO8uGUDFAr0zjqcmQuPH5Jm4Y7neHnrk8yaPINDOw6hPZukNdWMSH3LdyAp8NqWXykN9JcHZpf6Y9A2/X2U+vqG7IvnfXEs+xvog/50eEB/XwwdyM4pZcMFBiLdBBiXN6xSKWZXAcjdSaz3d9vz/LMUWtqYtKDyK3z01DJATQEqB5J0Ae0HmG7MzJw+h2RrhJkd3Xv2pZvfcQjzOw6pVXZE6lp0jx/7loIIWOm4t6wm1d8Ppdz2wEA6jVa6XaqynT5Syp73l6feiscYaxfHS3vsL0/TVSqlASZ3DqXyOXEse87gPkoTIuAWJ08f0+vXMkDtAiZV7GsHdh5gujEzo/0gTjrsbcxoP4iOtmm1elkR2YskKUBTgWQCtVGUgxeUg1h53F0pC2JZkIvn5cCWBT9KaXUrH/RKueuwx3lDBr5XjvPLrpk7t1QZTEslmqYePKEC1Ebgoop9i4FVw6Sbb2aT3P3VXLqNY5y/spZiM4tmLqzVy4lIg4r7TFntMx7VbWJQLQPUOiAxs4uAG4lefG8C7soncnc3s0eB5WZ2CfB24DTgbYiISMOo2Yq6aY+89xKBaTvRI+/97v6CmS01s3wT3hnAG4gxUDcD57n7E7XKq4iIjL+aDtRNg8wJw+y/Hbg9t72FCGYiItKgalaDEhER2R8KUCIiUpcUoEREpC5NxOmtiwCdnZ3jnQ8REaki9x1ddZT3RAxQcwGWLl063vkQEZF9mws8PdyBiRig1gPvALbCfi+RIyIitVEkgtMeSy5lkj2msBAREakD6iQhIiJ1SQFKRETqkgKUiIjUJQUoERGpSwpQIiJSlxSgRESkLilAiYhIXVKAEhGRujQRZ5I4YGZ2JLCCWOl3M3Cuu1cd5TzRmNkpwDXA64nFIr/i7t82sxZiFeQzidk5vubuV49fTmvHzDqAx4DL3f2WRiwLM5sLfAv4faAb+I67f74RywLAzI4DbgAMeAG4xt1vbqTyMLNjgR+4++x0e6/v3cz+CPgSMXPEA8DH3P35fb2OalCptIDvAb4PdADLgfvNbNq4ZqxGzGwecCfw18T7/zBwtZm9G7iK+DAeDhwDLDOzj45XXmtsBXBIbrsRy+IeYuqwOcBxxHs+mwYsCzMrEOVxg7tPJz4nN6Y/bid8eZhZYmbnA/cDLblDVd+7mS0BVgIfA2YAvwLWjOT1FKAGnQQ0u/v17t7r7muAJ4GzxjdbNbMAWO3ud7n7QFpz/BfgeGAZsNzdX3L3Z4GvAh8fr4zWipktA6YBj+d2N1RZmNlbgYXAp929292fIT4r62iwskgdBMwGEjNLgBLQB/TQGOVxFfBJ4ods3t7e+0eA+9z9Z+7eDVwCHG9mr9/XiylADVoCPFWxbxNwxDjkpebc/V/d/RPZtpkdTEy6+zBRLd+YSz7hy8XMDgOuAM7N7eug8criaCJAX2lmz5nZ08DpwKs0Xlng7tuIpqzvAr3ERKefI2qYjVAeK9z9aOAX2Y4RfC6W5I+5exewhRGUje5BDZoCdFXs6wLaxyEv48rMpgP3Aj8HNqS782UzocvFzIrAbcBn3b3TzLJDU9LHhikLIPuh8gBRk1oM/Ji49wKNVRZZE183cDbRJP52YC2wI00yocvD3f9nmN37+lwc8HerAtSgXcCkin3twM5xyMu4MbNFRBv7RmApg2WSL5uJXi6fB9zd11bs35U+NlJZ7AZecfcr0+1HzexmokkHGqssAD4AHO/uF6fbD5jZShq3PGDfn4sD/m5VE9+gjcRNvrzFDK22Tmhm9k6i1nQ3cGZ6z+EloJOhZTPRy+VDwJlmtsPMdhBNEd8kOs40WllsAtrTTkSZJqAR/78AmAe0VuzrI2qUjVgejOA7Ysh3q5m1A/MZQdmoBjVoHXHj8yKijfkMorv5XeOaqxoxs8OBHwCXuvvXKw7fClxhZo8R1fXPAn9b4yzWjLsvzm+b2SPA9Wk38500UFkA/0R8+V5rZn9BfNGcR9wo30xjlQVE77WrzexPgJuAo4A/Bs4HfkPjlUdmb98Rq4GfmdlJwIPA1cDD7v6f+7qoalApd+8B3ksEpu3ApcD73f2FvZ44cVwATCU+fDtzf18GLgeeIHo1rifa3leMX1bHVUOVRdrr6kTi/tNW4v7T37j7nTRYWQC4+5NEM9/HiftOq4G/cvd7aMDyyKn63t39caKz0QrgReB3gQ+O5KJaUVdEROqSalAiIlKXFKBERKQuKUCJiEhdUoASEZG6pAAlIiJ1SQFKRETqkgbqiuSY2S0MTlsznKuIWd7XAVPdvSZT2aTzA/4b8NHhBjia2YXE3IELxjgfNwLr3f27Y/k6IqAalEilzxAzM88llpUAODa376vAv6fPdw1z/lj5NPDoSEbfj7EvAF8wsxnjnA9pAKpBieS4+8vAywBmNjPd/YK7d1YkrdweM2bWRqyh865avWY17v68mf0E+BRw5ThnRyY4BSiR/ZTOKVZu4jOzErGy6iXEXHW/IBZpuxg4B3gFuMTdb03PnwpcSyyPXQJ+CnymylIGEJPX7nD3J3J5OAb4OjFf5HqiVpfPY7YS8pHpazwE/Km7bzKzHwHb3P2cXPrLgZPd/UQzOx34IvA6Ynqjb7r7V3KXXwusNLPl7t478pIT2T9q4hMZHdcAf0YsiT4f+CURmI4hvtC/bWbZujnfIQLZu4l57krAP5pZtR+MpxJz4AHlxSR/DDxKTFa6Crgod/xQYj2vvycWi3sXsa5TFmRuA04zs/wSCB8GbjezOcD3gevSPF4MLDezk3Np/5lYuvvoEZSLyAFTgBIZHd9w93Xu/ggxK/xO4HPu7sDXiPVwDjOzhUSN6Gx3X5/Wis4BFgDvqXLttxCTcGbOItZputDdN7n7SuCW3PEmosPEte7+jLs/RMw2vSQ9fjeQEIEPMzuKmAz2DuC3gWZgi7v/l7vfAZxMTAQKlCeQ3ZzmS2TMqIlPZHT8Ove8C3jW3bOZmLvTx1bg0PS551bqhVjAzYjgVmkOMQt05o3A4xXNa/9BGuDc/WkzW2tmf0nMHL2YaOrbmh7vMrO1RKC8g1gd9h/cfbuZvQTcTtToNgM/BFa5+/9W5GkbMLtaYYiMBtWgREZH5b2YgSrpmtK0vwe8Ofe3CPi7KucMEDWevMrtnuyJmb2RWGjwrcAGYm2eyyrS3wacmt4PO4sISrh7yd0/kuZvVXqNh8zsnIrzi0B/lfyKjArVoERq6ymiCW1y2hyImU0mAsaXic4MlTqBWbntx4gVf1vdfXe676jc8WXEgnAfyHaY2YcYGtR+Qqx7djEwDbgvTXcksMzd/xx4BLjKzFYT96huzZ0/kxr2ZJTGpAAlUkPu7mZ2L7DKzC4gVqtdTnSu2FTltA1EE11mDXAFcJOZfYkITucTAQei+W2xmb0DeA44jVj1dVsuHwNm9j0iQK1J7yuRXuOTaVPfbcQ9qePIBSczm040Va4/oEIQGSE18YnU3jKiK/rdxJf8dOAUd99RJf0Pid5+QHms1h8QQeJhognvulz6G4iedvcRwe0PgU8As83sd3LpVgNtpM176bW3ECvGnk50zLgzzefy3HknELWnh/fjPYvsN62oK1LnzKwdeBZ4j7v/chSveyrR5X2eu1e7Zzbced8DNrr7F0crLyLDUQ1KpM65excxhumC0biemR1mZh8k7nmt2M/gNJeozX1jNPIisjcKUCL/P1wHvMkq+qYfoHlEj8H/Jma02B+XAZe5+/Z9phR5jdTEJyIidUk1KBERqUsKUCIiUpcUoEREpC4pQImISF1SgBIRkbr0fxa4LksJ9hm/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_results(results.S, results.I, results.R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
